RT ANALYSIS

####THIS CHUNK WILL NEED TO BE ADJUSTED TO READ FROM GITHUB DATA###
#THIS CHUNK MAKES MORE SENSE TO MOVE AFTER FUNCTIONS ARE GENERATED
#IT IS LOCATED HERE FOR VISIBILITY ONLY, ONCE MOVED TO GITHUB, RELOCATE TO A MORE LOGICAL
#LOCATION IN THE CODE AS DESIRED
getwd()
## [1] "C:/Users/jeffb/Desktop/Life/personal-projects/COVID"
#read in full DSHS data with county and TSA information
full.dshs<-read.csv("combined-datasets/county.csv")#change path as needed
tsa.dshs<-read.csv("combined-datasets/tsa.csv")#change path as needed
rt.data.clean<-function(covid.data)#note: mydata will be either full.dshs or tsa.dshs data
  {

library(dplyr)
library(tidyr)
library(reshape2)
  #check the variable types of covid.data 
  #(comment out when running, this is for testing only)
    #str(covid.data)
  
  #keep variables from this list if they are in the data set
      #(county, date, population, TSA_Name, cases_daily, Population_DSHS, Population_Census)
  covid.data<-covid.data[,names(covid.data)%in% 
                           c("County", #name of texas county 
                             "Date", #date of data recorded
                             "TSA_Name", #name of texas TSA group for county
                             "TSA",
                             "Metro_Area", #indicateor of metro/nonmetro county
                             "Cases_Daily", #daily new cases
                             "Population_DSHS")] #population for county per Census estimates
                                                  #Note: Census data is from 2010, these values are projections 
                                                        #from that Census as provided by the Census Bureau
  #change Date variable to date format
  covid.data$Date<-as.Date(covid.data$Date)
  #change cases_daily to numeric - if not already numeric
  covid.data$Cases_Daily<-as.numeric(covid.data$Cases_Daily)
  #change any cases_daily below zero to zero
  #NOTE: We shouldn'te actually have cases_daily below zero, 
   #these are a result of data problems related to how cases are counted
   #we will have to address this in our data as this will cause 
   #potential issues with rt estimates
  covid.data<-covid.data%>%mutate(Cases_Daily=ifelse(Cases_Daily<0, 0, Cases_Daily))
  return(covid.data)
}
#separate county, metro, TSA and State data into separate Data frames
rt.data.organize<-function(mycounty, mytsa){
  library(dplyr)
  library(tidyr)
  #call the rt.data.clean function on mydata
  cleaned.county<-rt.data.clean(mycounty)
  cleaned.tsa<-rt.data.clean(mytsa)

  #create metro.df by mutating by grouping by Metro
  metro.temp<-cleaned.county%>%
    group_by(Metro_Area,Date)%>%
    mutate(metro_cases_daily=sum(Cases_Daily, na.rm=TRUE))%>%#generate variable of state daily cases
    mutate(metro_pop_DSHS=sum(Population_DSHS, na.rm=TRUE))%>%#generate state population_DSHS variable
    dplyr::select(Date,Metro_Area, metro_cases_daily, metro_pop_DSHS)%>%#select only variables of interest
    distinct()#keep only distinct rows (this will drop repeated rows, now that we have only state level va
  metro.df<-data.frame(Date=metro.temp$Date,
                       Metro_Area=metro.temp$Metro_Area,
                       Cases_Daily=metro.temp$metro_cases_daily,
                       Population_DSHS=metro.temp$metro_pop_DSHS)
  #Drop TSA_Name from county data frame
  county.df<-cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area)
  
   #create TSA.df, use rt.data.clean function
  TSA.df<-rt.data.clean(mytsa)
  
  #create state.df
  #at the state level we want daily new cases, and population variables from DSHS and Census data
  state.temp<-cleaned.tsa%>%
    group_by(Date)%>%
    mutate(state_cases_daily=sum(Cases_Daily, na.rm=TRUE))%>%#generate variable of state daily cases
    mutate(state_pop_DSHS=sum(Population_DSHS, na.rm=TRUE))%>%#generate state population_DSHS variable
    dplyr::select(Date,state_cases_daily, state_pop_DSHS)%>%#select only variables of interest
    distinct()#keep only distinct rows (this will drop repeated rows, now that we have only state level vars)
  state.df<-data.frame(Date=state.temp$Date, 
                       Cases_Daily=state.temp$state_cases_daily,
                       Population_DSHS=state.temp$state_pop_DSHS)
  #return a list of the three dataframes generated
 return(list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df))
}
rt.df.extraction<-function(Rt.estimate.output){

   #first convert named vector back to dataframe
  #reorder and rename columns from the R0.estimate output
  rt.df<-setNames(stack(Rt.estimate.output$estimates$TD$R)[2:1], c('Date', 'Rt'))
 
  #extract the upper and lower limits of the Rt 95% confidence intervals, these will be lists
  CI.lower.list<-Rt.estimate.output$estimates$TD$conf.int$lower
  CI.upper.list<-Rt.estimate.output$estimates$TD$conf.int$upper

  #use unlist function to format as vector
  CI.lower <- unlist(CI.lower.list, recursive = TRUE, use.names = TRUE)
  CI.upper <- unlist(CI.upper.list, recursive = TRUE, use.names = TRUE)

  #add the upper and lower bounds to the dataframe
  rt.df$lower<-CI.lower
  rt.df$upper<-CI.upper
  
  rt.df$Date = as.Date(rt.df$Date)
  
  return(rt.df)
}
#############################
#Remove this code later, for testing purposes only
#creating test county
#testcounty.df<-county.daily%>%subset(County=="Archer")
#testcounty.rt<-covid.rt(testcounty.df)
  #mydata<-testcounty.df

#Test the covid.function on TSA subset
#  #TSA.G<-TSA.daily%>%subset(TSA_Name=="Panhandle RAC")

###############################
#Function to generate Rt estimates, need to read in data frame and population size
covid.rt<-function(mydata){
library(R0)
  #change na values to 0
  mydata<-mydata%>%replace(is.na(.),0)
  #create a variable that is the sum of Cases_Daily
  sum.daily.cases<-sum(mydata$Cases_Daily)
  
  #create a vector of new cases
  mydata.new<-pull(mydata, Cases_Daily)
  #create a vector of dates
  date.vector<-pull(mydata, Date)
  #create a named numerical vector using the date.vector as names of new cases
    #Note: this is needed to run R0 package function estimate.R()
  names(mydata.new)<-c(date.vector)
  
  #create the generation time 
  ##### NOTE: THIS MAY NEED MODIFICATION ONCE MORE INFORMATIONIS AVAILABLE ####

  #gen.time<-generation.time("lognormal", c(4.0, 2.9))#verify these values and distribution/update as needed
  #gen.time<-generation.time("lognormal", c(4.7,2.9)) #Nishiura

  #Tapiwa, Ganyani "Esimating the gen interval for Covid-19":
  # gen.time<-generation.time("gamma", c(5.2, 1.72)) #Singapore
  #gen.time<-generation.time("gamma", c(3.95, 1.51)) #Tianjin
  gen.time<-generation.time("gamma", c(3.96, 4.75))

  #get DSHS population and census population
  pop.DSHS<-mydata$Population_DSHS[1]

  #Find min.date, max.date and row number of max.date
  min.date<-min(mydata$Date)
  max.date<-max(mydata$Date)
  # str(mydata)

  #get row number of March 9 and first nonzero entry
  #find max row between the two (this will be beginning of rt data used)
  march9.row<-which(mydata$Date=="2020-03-09")
  first.nonzero<-min(which(mydata$Cases_Daily>0))
  last.nonzero<-max(which(mydata$Cases_Daily>0))
  minrow<-max(march9.row,first.nonzero)
  
  #Set last row that will be used, either todays date or last nonzero row
  j<-as.integer(min(last.nonzero, max.date))
  
  if(is.na(minrow)){
    rt.DSHS.df<-data.frame(Date=as.Date(mydata$Date), Rt=rep(NA, j), lower=rep(NA, j), upper=rep(NA, j))
  }
  
  if(sum.daily.cases<=50){
    k=length(mydata$Date)
    rt.DSHS.df<-data.frame(Date=as.Date(mydata$Date), Rt=rep(NA, k), lower=rep(NA, k), upper=rep(NA, k))
  }
  #only move forward if both previous error catch options are false
  if(is.na(minrow)==FALSE & sum.daily.cases>50){
    #get the row number for the mindrow date
    i<-minrow
  
    #reduce the vector to be rows from min date (March 9 or first nonzero case) to current date
    mydata.newest<-mydata.new[i:j]
    
    
    rt.DSHS<-estimate.R(mydata.newest, 
                        gen.time, 
                        #begin=as.integer(1), 
                        #end=total.days, 
                        methods=c("TD"), 
                        pop.size=pop.DSHS,
                        nsim=1000)
  
      rt.DSHS.df<-rt.df.extraction(rt.DSHS)
  }
  
  rt.out = subset(rt.DSHS.df, Date > as.Date('2020-03-09'))
  return(rt.out)
}

County

Metro

TSA

#run covid.rt function on state level data 
#(this will be a data frame output)
state.rt.df<-covid.rt(state.daily)
## Warning in est.R0.TD(epid = c(`2020-03-09` = 7, `2020-03-10` = 3, `2020-03-11` =
## 3, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-09` = 7, `2020-03-10` = 3, `2020-03-11` =
## 3, : Using initial incidence as initial number of cases.
########## KABLE TABLES CAN BE COMMENTED OUT WHEN RUN IN GIT ##############
# library(kableExtra)
# county.rt.df[1:25,]%>%kable(caption="Texas County Level RT Estimates")%>%kable_styling(full_width=FALSE)
# metro.rt.df[1:25,]%>%kable(caption="Texas Metro Level RT Estimates")%>%kable_styling(full_width=FALSE)
# TSA.rt.df[1:25,]%>%kable(caption="Texas TSA RT Estimates")%>%kable_styling(full_width = FALSE)
# state.rt.df[1:25,]%>%kable(caption="Texas State RT Estimates")%>%kable_styling(full_width=FALSE)

#NOTE: FILE DESTINATION WILL NEED TO BE UPDATED FOR EACH OF THESE TO 
  #MATCH DIRECTORY IN YOUR WORKING DIRECTORY

write.csv(county.rt.df, file="statistical-output/rt/county_rt.csv", row.names=FALSE)
write.csv(metro.rt.df, file="statistical-output/rt/metro_rt.csv", row.names=FALSE)
write.csv(TSA.rt.df, file="statistical-output/rt/tsa_rt.csv", row.names=FALSE)
write.csv(state.rt.df, file="statistical-output/rt/state_rt.csv", row.names=FALSE)

TIME SERIES

#ts.rt.data.clean function will clean data, dropping unwanted variables
#note: mydata will be either full.dshs or tsa.dshs data
ts.data.clean<-function(covid.data) {
  
  library(dplyr)
  library(tidyr)
  library(reshape2)
  # str(covid.data)
  #keep variables from this list if they are in the data set
      #(county, date, population, TSA_Name, Cases_Daily, Population_DSHS, Population_Census)
  covid.data<-covid.data[,names(covid.data)%in% 
                           c("County", #name of texas county 
                             "Date", #date of data recorded
                             "TSA",
                             "TSA_Name", #name of texas TSA group for county
                             "Metro_Area", #indicateor of metro/nonmetro county
                             "Cases_Daily", #daily new cases
                             "Cases_Cumulative",
                             "Population_DSHS")] #population for county per DSHS data
                                                  #Note: Census data is from 2010, these values are projections 
                                                        #from that Census as provided by the Census Bureau
  #change Date variable to date format
  covid.data$Date<-as.Date(covid.data$Date)
  #change Cases_Daily to numeric - if not already numeric
  covid.data$Cases_Daily<-as.numeric(covid.data$Cases_Daily)
  #change any Cases_Daily below zero to zero
  covid.data<-covid.data%>%mutate(Cases_Daily=ifelse(Cases_Daily<0, 0, Cases_Daily))
  return(covid.data)
  
}
#separate county, metro, TSA and State data into separate Data frames
ts.data.organize<-function(mycounty, mytsa){
  library(dplyr)
  library(tidyr)
  
  ### COUNTY ###
  
  #call the ts.data.clean function on mydata
  cleaned.county<-ts.data.clean(mycounty)
  
  # drop NA dates
  cleaned.county<-subset(cleaned.county, !is.na(Date) & Date>=as.Date('2020-03-04'))

  #Drop TSA_Name from county data frame
  county.df<-cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area)
  county.df$Level = 'County'
  colnames(county.df)[1] = 'Level_Name'
  ### METRO ###

  #create metro.df by mutating by grouping by Metro
  metro.temp<-cleaned.county%>%
    group_by(Metro_Area,Date)%>%
    mutate(metro_Cases_Daily=sum(Cases_Daily, na.rm=TRUE))%>%#generate variable of state daily cases
    mutate(metro_Cases_Cumulative=sum(Cases_Cumulative, na.rm=TRUE))%>%#generate variable of state daily cases
    mutate(metro_pop_DSHS=sum(Population_DSHS, na.rm=TRUE))%>%#generate state population_DSHS variable
    dplyr::select(Date,Metro_Area, metro_Cases_Daily, metro_Cases_Cumulative, metro_pop_DSHS)%>%#select only variables of interest
    distinct()#keep only distinct rows (this will drop repeated rows, now that we have only state level va
    
  metro.df<-data.frame(Date=metro.temp$Date,
                       Level = 'metro',
                       Level_Name=metro.temp$Metro_Area,
                       Cases_Daily=metro.temp$metro_Cases_Daily,
                       Cases_Cumulative = metro.temp$metro_Cases_Cumulative,
                       Population_DSHS=metro.temp$metro_pop_DSHS)
  
    
  # drop NA dates
  metro.df<-subset(metro.df, !is.na(Date) & Date>=as.Date('2020-03-04'))
  
  ### TSA ###
  TSA.df<-ts.data.clean(mytsa)
  TSA.df<-subset(TSA.df, Date>=as.Date('2020-03-04'))
  
  TSA.df$Level = 'TSA'
  colnames(TSA.df)[2] = 'Level_Name'
  ### STATE ###
  
  #create state.df
  #at the state level we want daily new cases, and population variables from DSHS and Census data
  state.temp<-TSA.df%>%
    group_by(Date)%>%
    mutate(state_Cases_Daily=sum(Cases_Daily, na.rm=TRUE))%>%#generate variable of state daily cases
    mutate(state_Cases_Cumulative=sum(Cases_Cumulative, na.rm=TRUE))%>%#generate variable of state Cumulative cases
    mutate(state_pop_DSHS=sum(Population_DSHS, na.rm=TRUE))%>%#generate state population_DSHS variable
    dplyr::select(Date,state_Cases_Daily, state_Cases_Cumulative, state_pop_DSHS)%>%#select only variables of interest
    distinct()#keep only distinct rows (this will drop repeated rows, now that we have only state level vars)
    
  state.df<-data.frame(Date=state.temp$Date,
                       Level = 'State',
                       Level_Name = 'Texas',
                       Cases_Daily=state.temp$state_Cases_Daily,
                       Cases_Cumulative=state.temp$state_Cases_Cumulative, 
                       Population_DSHS=state.temp$state_pop_DSHS)
  
  state.df<- subset(state.df, Date>=as.Date('2020-03-04'))
  
  #return a list of the three dataframes generated
  return(list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df))
}
library(forecast)
## Warning: package 'forecast' was built under R version 3.6.3
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:nlme':
## 
##     getResponse
library(ggplot2) 

library(astsa)
## Warning: package 'astsa' was built under R version 3.6.3
## 
## Attaching package: 'astsa'
## The following object is masked from 'package:forecast':
## 
##     gas
library(forecast)
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(ggplot2)
library(fpp2)
## Warning: package 'fpp2' was built under R version 3.6.3
## Loading required package: fma
## Warning: package 'fma' was built under R version 3.6.3
## 
## Attaching package: 'fma'
## The following objects are masked from 'package:astsa':
## 
##     chicken, sales
## The following objects are masked from 'package:MASS':
## 
##     cement, housing, petrol
## Loading required package: expsmooth
## Warning: package 'expsmooth' was built under R version 3.6.3
## 
## Attaching package: 'fpp2'
## The following object is masked from 'package:astsa':
## 
##     oil
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
covid.arima.forecast<-function(mydata)
{
  maxdate<-max(mydata$Date)
  mindate <-min(mydata$Date)
  prediction.period<-10 #this can be changed as desired
  
  if(max(mydata$Cases_Daily>=5, na.rm = TRUE))
  {
    
    #for the time series, you need a start day, which is the year, and the day number - 
    #we need to double check how this works exactly
    my.timeseries<-ts(mydata$Cases_Daily)
    
    #d=0 restricts first differencing to 0 so that daily cases aren't differenced
    #get a fit model using the best ARIMA from auto.arima
    arima.fit<-forecast::auto.arima(my.timeseries, d=0)
    
    #use predict() on the fit model
    # prediction.period is how many more days to forecast
    # 10 day forecast
    arima.forecast<-forecast::forecast(arima.fit, h = prediction.period) 
    
    #auto correlation function
    # diagnostics
    acf<-ggAcf(my.timeseries, lag.max=30)
    pacf<-ggPacf(my.timeseries, lag.max=30)
    grid.arrange(acf, pacf, nrow=2)
    ggsave(paste0('statistical-output/time-series/diagnostics/', mydata$Level[1],'/', mydata$Level_Name[1], '.png'))
     
    #return a dataframe of the arima model(Daily cases by date) 
    #and forecasted values ***NOT SURE HOW THIS WILL WORK***
    arima.out<-data.frame(Date = seq(mindate, maxdate + 10, by = 'days'),
                          Cases_Daily = c(my.timeseries, arima.forecast[['mean']]),
                          CI_Lower = c(rep(NA, times = length(my.timeseries)), arima.forecast[['lower']][, 2]),
                          CI_Upper = c(rep(NA, times = length(my.timeseries)), arima.forecast[['upper']][, 2]))

    } else {
      
    arima.out<-data.frame(Date = seq(mindate, maxdate + 10, by = 'days'),
                          Cases_Daily = c(mydata$Cases_Daily, rep(NA, times = 10)),
                          CI_Lower = rep(NA, times = length(mydata$Date) + 10),
                          CI_Upper =  rep(NA, times = length(mydata$Date) + 10)
                          )
  }
  return(arima.out)
}

County

#get county forecasts using covid.arima.forecast function

#load nlme package
library(nlme)
#apply the covid.arima.forecast function to the county.Daily dataframe by county
#note that the output will be a list of dataframes
county.arima.output<-nlme::gapply(county.Daily, FUN = covid.arima.forecast, groups=county.Daily$Level_Name)
## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

#load data.table package
library(data.table)
county.arima.df<-data.table::rbindlist(county.arima.output, idcol=TRUE)

colnames(county.arima.df)[1] = 'County'

Metro

#get metro forecasts using covid.arima.forecast function
#load nlme package
library(nlme)
#apply the covid.arima.forecast function to the metro.Daily dataframe by metro
#note that the output will be a list of dataframes
metro.arima.output<-nlme::gapply(metro.Daily, FUN = covid.arima.forecast, groups=metro.Daily$Level_Name)
## Saving 7 x 5 in image

## Saving 7 x 5 in image

#load data.table package
library(data.table)
metro.arima.df<-data.table::rbindlist(metro.arima.output, idcol=TRUE)

colnames(metro.arima.df)[1] = 'Metro_Area'

TSA

#ge#load nlme package
library(nlme)
#apply the covid.arima.forecast function to the TSA.Daily dataframe by TSA
#note that the output will be a list of dataframes
TSA.arima.output<-nlme::gapply(TSA.Daily, FUN = covid.arima.forecast, groups=TSA.Daily$Level_Name)
## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

#load data.table package
library(data.table)
TSA.arima.df<-data.table::rbindlist(TSA.arima.output, idcol=TRUE)


colnames(TSA.arima.df)[1] = 'TSA'
TSA.arima.df = merge(TSA.arima.df, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = c('TSA'))
#get Texas forecasts using covid.arima.forecast function
texas.arima.df<-covid.arima.forecast(state.Daily)
## Saving 7 x 5 in image

#export arima results into .csv files
write.csv(texas.arima.df, 'statistical-output/time-series/state_timeseries.csv', row.names = F)
write.csv(TSA.arima.df, 'statistical-output/time-series/tsa_timeseries.csv', row.names = F)
write.csv(county.arima.df, 'statistical-output/time-series/county_timeseries.csv', row.names = F)
write.csv(metro.arima.df, 'statistical-output/time-series/metro_timeseries.csv', row.names = F)

STANDARD STATISTICAL TESTS

Ratio of cases

TSA level

library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following object is masked from 'package:base':
## 
##     date
cumcases.TSA <- TSA.Daily %>% dplyr::select(Date,Level_Name, Cases_Cumulative)
cumcases.TSA$Date <- ymd(cumcases.TSA$Date)

latestdate <- max(cumcases.TSA$Date)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.TSA, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.TSA, Date==middate, select=Cases_Cumulative)

week1 <- subset(cumcases.TSA, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.TSA, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
#current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
tsa_current_ratio_dat <- data.frame(TSA=subset(cumcases.TSA, Date==latestdate, select=Level_Name)[,1],
                                    current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

# add TSA name
tsa_current_ratio_dat = merge(tsa_current_ratio_dat, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = 'TSA')

County level

cumcases.county <- county.Daily %>% dplyr::select(Date,Level_Name, Cases_Cumulative)
cumcases.county$Date <- ymd(cumcases.county$Date)

latestdate <- max(cumcases.county$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.county, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.county, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.county, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.county, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
county_current_ratio_dat <- data.frame(County=subset(cumcases.county, Date==latestdate, select=Level_Name)[,1], current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

Metropolitan level

cumcases.metro <- metro.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative)
cumcases.metro$Date <- ymd(cumcases.metro$Date)

latestdate <- max(cumcases.metro$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.metro, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.metro, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.metro, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.metro, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA
metro_current_ratio_dat <- data.frame(Metro_Area=subset(cumcases.metro, Date==latestdate, select=Level_Name)[,1], current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

Texas level

cumcases.state <- state.Daily %>% dplyr::select(Date, Cases_Cumulative)
cumcases.state$Date <- ymd(cumcases.state$Date)

latestdate <- max(cumcases.state$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.state, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.state, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.state, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.state, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA
state_current_ratio_dat <- data.frame(current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

export

#export arima results into .csv files
write.csv(state_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/state_case_ratio.csv', row.names = F)
write.csv(tsa_current_ratio_dat,'statistical-output/standard-stats/case-ratios/tsa_case_ratio.csv', row.names = F)
write.csv(county_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/county_case_ratio.csv', row.names = F)
write.csv(metro_current_ratio_dat,'statistical-output/standard-stats/case-ratios/metro_case_ratio.csv', row.names = F)

rolling average

Omitted fow now, Dr. Yaseen is processing this in Tableau

MACHINE LEARNING

TODO: save output

County

set.seed(1)
train= sample(1:dim(datatrend.dat)[1],dim(datatrend.dat)[1]/2)
test=-train
datatrend.county.train=datatrend.dat[train,]    #traiing set
datatrend.county.test=datatrend.dat[test,]     #testing set

#logistic
glm.county.fit=glm(ratio2 ~ MKteststat+ MKpvalue+ Spearteststat+ Spearpvalue +Population + black_per +hispanic_per + age0_4 + age10_14 +age15_19+age20_24+  age25_29 + age30_34 + age35_39 + age40_44 + age45_49  +  age5_9+  age50_54+  age55_59 + age60_64 + age65_69+ age70_74+ age75_79+ age80_84 +  age85up, data=datatrend.county.train,family=binomial)     #fit logistic model
summary(glm.county.fit)    #summany of logistics model
## 
## Call:
## glm(formula = ratio2 ~ MKteststat + MKpvalue + Spearteststat + 
##     Spearpvalue + Population + black_per + hispanic_per + age0_4 + 
##     age10_14 + age15_19 + age20_24 + age25_29 + age30_34 + age35_39 + 
##     age40_44 + age45_49 + age5_9 + age50_54 + age55_59 + age60_64 + 
##     age65_69 + age70_74 + age75_79 + age80_84 + age85up, family = binomial, 
##     data = datatrend.county.train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6960  -1.2005   0.9214   1.1302   1.7274  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.908e-01  1.173e-01   3.331 0.000867 ***
## MKteststat     1.540e+00  1.364e+00   1.129 0.258796    
## MKpvalue      -9.409e-01  5.931e-01  -1.586 0.112647    
## Spearteststat -5.468e-01  1.090e+00  -0.502 0.615905    
## Spearpvalue    8.108e-01  6.168e-01   1.315 0.188676    
## Population     3.267e-07  6.660e-07   0.491 0.623773    
## black_per     -6.244e-03  6.691e-03  -0.933 0.350655    
## hispanic_per  -2.926e-03  1.813e-03  -1.614 0.106615    
## age0_4        -1.416e-04  1.332e-04  -1.063 0.287692    
## age10_14      -2.969e-04  2.221e-04  -1.336 0.181406    
## age15_19       2.981e-04  1.303e-04   2.289 0.022106 *  
## age20_24      -8.398e-05  5.175e-05  -1.623 0.104609    
## age25_29       8.228e-05  1.234e-04   0.667 0.505042    
## age30_34      -2.097e-04  2.320e-04  -0.904 0.365961    
## age35_39       4.172e-04  2.760e-04   1.512 0.130601    
## age40_44      -3.080e-04  2.502e-04  -1.231 0.218324    
## age45_49      -1.861e-04  2.709e-04  -0.687 0.492167    
## age5_9         2.379e-04  2.595e-04   0.916 0.359437    
## age50_54       3.448e-04  3.194e-04   1.080 0.280313    
## age55_59      -9.687e-05  2.842e-04  -0.341 0.733210    
## age60_64      -8.505e-05  2.233e-04  -0.381 0.703287    
## age65_69       1.873e-04  2.767e-04   0.677 0.498324    
## age70_74      -1.975e-04  4.094e-04  -0.482 0.629554    
## age75_79       1.724e-05  5.619e-04   0.031 0.975519    
## age80_84       1.400e-05  5.961e-04   0.023 0.981262    
## age85up       -1.480e-04  3.235e-04  -0.457 0.647367    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4793.4  on 3471  degrees of freedom
## Residual deviance: 4726.8  on 3446  degrees of freedom
##   (516 observations deleted due to missingness)
## AIC: 4778.8
## 
## Number of Fisher Scoring iterations: 4
glm.county.probs=predict(glm.county.fit, datatrend.county.test, type="response")     #predict
glm.county.predict=rep("0",3399)
glm.county.predict[glm.county.probs > 0.5]="1"      #round-off
#glm.county.predict=as.factor(glm.county.predict)
table(glm.county.predict, datatrend.county.test$ratio2)
##                   
## glm.county.predict    0    1
##                  0  428  924
##                  1 1073 1319
mean(glm.county.predict==datatrend.county.test$ratio2)
## [1] NA
#random forest
#rf.county.train=randomForest(ratio ~ MKteststat+ MKpvalue+ Spearteststat+ Spearpvalue +Population + black_per +hispanic_per + age0_4 + age10_14 +age15_19+age20_24+  age25_29 + age30_34 + age35_39 + age40_44 + age45_49  +  age5_9+  age50_54+  age55_59 + age60_64 + age65_69+ age70_74+ age75_79+ age80_84 +  age85up, data=datatrend.county.train, mtry=5, importance =TRUE)
#yhat.rf.county = predict(rf.county.train, newdata=datatrend.dat.rf[-train ,])
#mean((yhat.rf.county-datatrend.county.test)^2)

TSA

set.seed(1)
train= sample(1:dim(tsa_datatrend.dat)[1],dim(tsa_datatrend.dat)[1]/2)
test=-train
datatrend.tsa.train=tsa_datatrend.dat[train,]    #traiing set
datatrend.tsa.test=tsa_datatrend.dat[test,]     #testing set

#logistic
glm.tsa.fit=glm(ratio2 ~ MKteststat+ MKpvalue+ Spearteststat+ Spearpvalue +Population + black_per +hispanic_per + Retail_Recreation+ Grocery_Pharmacy   +  Parks  +  Transit +Workplaces+ Residential+age0_4 + age10_14 +age15_19+age20_24+  age25_29 + age30_34 + age35_39 + age40_44 + age45_49  +  age5_9+  age50_54+  age55_59 + age60_64 + age65_69+ age70_74+ age75_79+ age80_84 +  age85up, data=datatrend.tsa.train ,family=binomial)     #fit logistic model
summary(glm.tsa.fit)    #summany of logistics model
## 
## Call:
## glm(formula = ratio2 ~ MKteststat + MKpvalue + Spearteststat + 
##     Spearpvalue + Population + black_per + hispanic_per + Retail_Recreation + 
##     Grocery_Pharmacy + Parks + Transit + Workplaces + Residential + 
##     age0_4 + age10_14 + age15_19 + age20_24 + age25_29 + age30_34 + 
##     age35_39 + age40_44 + age45_49 + age5_9 + age50_54 + age55_59 + 
##     age60_64 + age65_69 + age70_74 + age75_79 + age80_84 + age85up, 
##     family = binomial, data = datatrend.tsa.train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4102  -0.8869   0.4935   0.8166   1.6311  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        1.260e+01  1.538e+01   0.819  0.41283   
## MKteststat        -6.272e+00  4.306e+00  -1.457  0.14518   
## MKpvalue          -2.144e+00  1.635e+00  -1.312  0.18968   
## Spearteststat      5.578e+00  3.263e+00   1.709  0.08738 . 
## Spearpvalue        1.959e+00  1.654e+00   1.185  0.23618   
## Population         1.952e-05  5.719e-05   0.341  0.73285   
## black_per          3.232e-01  2.264e+00   0.143  0.88646   
## hispanic_per      -1.862e-01  3.074e-01  -0.606  0.54464   
## Retail_Recreation  1.014e-01  7.921e-02   1.280  0.20045   
## Grocery_Pharmacy  -5.747e-02  5.368e-02  -1.071  0.28437   
## Parks             -1.261e-02  8.397e-03  -1.502  0.13313   
## Transit           -6.054e-02  2.204e-02  -2.747  0.00602 **
## Workplaces        -1.407e-01  9.186e-02  -1.532  0.12557   
## Residential       -2.589e-01  3.011e-01  -0.860  0.38976   
## age0_4             3.629e-03  1.460e-02   0.249  0.80364   
## age10_14           5.805e-03  1.994e-02   0.291  0.77089   
## age15_19          -5.806e-03  1.210e-02  -0.480  0.63144   
## age20_24           2.540e-03  5.831e-03   0.436  0.66312   
## age25_29          -3.754e-03  1.359e-02  -0.276  0.78235   
## age30_34           1.909e-03  1.174e-02   0.163  0.87076   
## age35_39           6.846e-04  6.402e-03   0.107  0.91484   
## age40_44          -9.738e-03  1.414e-02  -0.689  0.49090   
## age45_49           1.887e-02  3.942e-02   0.479  0.63210   
## age5_9            -4.514e-03  2.255e-02  -0.200  0.84137   
## age50_54          -1.155e-02  3.032e-02  -0.381  0.70338   
## age55_59           6.267e-04  4.141e-03   0.151  0.87971   
## age60_64           4.545e-04  3.774e-03   0.120  0.90414   
## age65_69          -1.490e-03  6.521e-03  -0.228  0.81928   
## age70_74           7.705e-03  2.070e-02   0.372  0.70972   
## age75_79          -2.510e-02  6.635e-02  -0.378  0.70521   
## age80_84           2.725e-02  6.913e-02   0.394  0.69348   
## age85up           -3.031e-03  6.396e-03  -0.474  0.63558   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 330.58  on 270  degrees of freedom
## Residual deviance: 277.47  on 239  degrees of freedom
##   (196 observations deleted due to missingness)
## AIC: 341.47
## 
## Number of Fisher Scoring iterations: 5
glm.tsa.probs=predict(glm.tsa.fit,datatrend.tsa.test,type="response")     #predict
glm.tsa.predict=rep("0",401)
glm.tsa.predict[glm.tsa.probs > 0.5]="1"      #round-off
#glm.county.predict=as.factor(glm.county.predict)
table(glm.tsa.predict, datatrend.tsa.test$ratio2)
##                
## glm.tsa.predict   0   1
##               0  84 141
##               1  59 153
mean(glm.tsa.predict==datatrend.tsa.test$ratio2)
## [1] NA